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Abstract 

A semi-spectral Chebyshev method for solving numerically singular integral equa- 
tions is presented and applied in the quarkonium bound-state problem in momentum 
space. The integrals containing both, logarithmic and Cauchy singular kernels, can 
be evaluated without subtractions by dedicated automatic quadratures. By intro- 
ducing a Chebyshev mesh and using the Nystrom algorithm the singular integral 
equation is converted into an algebraic eigenvalue problem that can be solved by 
standard methods. The proposed scheme is very simple to use, is easy in program- 
ming and highly accurate. 
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1 Introduction 

In a recent work [1] we have advocated the Chebyshev semi-spectral method 
demonstrating its efficiency in solving some typical differential and integral 
equations emerging in quantum mechanics. The present paper is in the same 
vein but here we wish to focus our attention solely on the heavy quarkonium 
momentum space bound-state problem. Admittedly, the problem is not new 
but our incentive here is to examine the effectiveness of the semi-spectral 
approach in solving strongly singular integral equations. Since the latter topic 
was beyond the scope of [1], this work may be regarded as an immediate 
continuation of the previous paper. 

We would like to believe that the presented method will be useful also out- 
side quantum mechanics, especially that strongly singular integral equations 
are encountered in many areas of science and engineering. The well known 
physical applications comprise the quantum mechanical scattering problem, 
the Omnes formulation [2] of the final-state-interaction, radiative transfer, 
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neutron transport [3] etc. The list of engineering applications is by no means 
restricted to the widely known aerofoil problem [4] and, indeed, many im- 
portant problems of engineering mechanics like elasticity, plasticity, fracture 
mechanics, etc. may be also efficiently expressed in terms of singular and hy- 
persingular integral equations. Because it is not always possible to find explicit 
solutions to the problems posed, much attention has been devoted to approx- 
imate methods. It is interesting to note that even when an analytic solution 
is known, quite often the latter takes the form of a singular integral whose 
numerical evaluation might be more complicated than a numerical solution of 
the integral equation. 

A hypersingular integral equation arises in quantum mechanics already at a 
quite elementary level when the linear potential bound-state problem, easily 
tackled in configuration space, is approached in momentum space. This prob- 
lem is far from academic since the linear potential plays an important role 
not only in atomic physics where it is associated with the hydrogen radial 
Stark effect but also in particle physics serving as a simple confinement model 
of QCD. Although, in principle, QCD alone should describe the spectroscopy 
of heavy quarkonia but the implementation of such program is very difficult 
and instead various phenomenological models incorporating some QCD prop- 
erties have been developed (for a recent review of quarkonium physics and 
references to the literature cf. [5]). The QCD motivated quark potential mod- 
els have played a prominent role in understanding quarkonium spectroscopy 
and are capable of reproducing with surprising accuracy a sizable part of the 
meson and baryon properties. The non-relativistic potential approach may be 
justified by the fact that the bottom quark and, perhaps to a lesser extent, 
also the charmed quark have masses that are large in comparison with A - the 
typical QCD hadronic mass scale. The quark-ant iquark potential has been 
tailored to mock up the properties expected from QCD and the different po- 
tential shapes set up in the early days after years of research have evolved to 
a common form that one might expect from the asymptotic limits of QCD. 
The prototype for these potentials is still the popular Cornell potential [6] 
including the one-gluon-exchange Coulomb potential supplemented by a lin- 
ear potential simulating confinement, as expected from QCD. Therefore, this 
potential will be also considered in this paper. 

Obviously, the non-relativistic potential model can not be pushed beyond cer- 
tain limits and for systems containing one light quark a complete disregard of 
relativistic effects might be a serious omission. In addition to that, it was some- 
what embarrassing when people realized [7] that within the non-relativistic 
formalism the mesons containing a light quark might be more massive than 
a meson composed with heavier quarks. These difficulties could be aleviated 
at the expense of a semirelativistic treatment where the relativistic expres- 
sion for the energy is used. A popular relativistic extension of the Schrodinger 
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equation is the spinless Salpeter equation 



\jp 2 + m\ + \jp 2 + m 2 2 + V(r) *(r) = E *(r) 



(1) 



where mi, m 2 are the quark masses, p is the cm. momentum, V(r) denotes 
the quark- ant iquark potential and E is the eigenenergy. Since in such case 
the Laplacian operator appears under a square root, the coordinate space 
is rather unwieldy for solving the bound state problem and the momentum 
space seems to be the most natural alternative. Indeed, in momentum space 
the energy operator is diagonal and the difference in computational effort 
between non-relativistic and semi-relativistic treatment is minor. Although 
the momentum space approach solves some problems automatically but at 
the same time it does create another difficulty in that the quark-antiquark 
potential gives rise to a singular kernel in the appropriate integral equation. 
Whilst the Coulomb potential yields a kernel with a logarithmic singularity 
that can be removed by subtraction [8], the kernel associated with a linear 
potential exhibits a double-pole singularity for which the subtraction scheme 
is insufficient. To clarify this important point let us consider just the linear 
potential for simplicity restricting our attention to a zero orbital momentum 
state. The potential term that enters the appropriate wave equation involves 
the integral with a double pole singularity 



where <f>(k) is the wave function, <f>'(p) denotes the derivative and p is a real 
parameter. It may be easily verified that the two extra terms occurring on 
the right hand side of (2) can be supplemented with impunity because the 
integrals multiplying, respectively, <j)(p) and <j)'(p) are both bound to vanish. 
The integral on the right hand side is non-singular and in the limit k — > p the 
integrand goes to a finite limit \4>'(p)/p + \4>"(p). This demonstrates explic- 
itly that by using a subtraction technique it is perfectly possible to remove 
the singularity converting the integral to a form amenable for computation. 
Nevertheless, the subtraction scheme (2) would be insufficient for solving an 
integral equation as it introduces unknown first <f>'(p) and second derivative 
4>"(p) at the top of the unknown function. This also explains why the Nys- 
trom method, which has been rather efficient in solving the Coulomb bound 
state problem in momentum space [8], does not work for the linear potential. 
( The calculation using Nystrom method presented in [9] is incorrect because 
the infinite diagonal term in the potential matrix has been simply omitted 
whereas the proposed correction, given in their eq. (34), is proportional to a 
logarithmically diverging integral.) 

In the early attempts to overcome this difficulty the singularity was removed 
by hand, by introducing an arbitrary cut-off [10] [11] in the potential. The re- 
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suiting non-singular integral equation involving the modified potential could 
be then solved by standard methods. The unwelcome arifacts of the cutoff 
might be eventually disposed of by perturbative methods [11]. However, a 
more promising approach is to seek the wave function in the form of an ex- 
pansion in terms of a complete set of orthogonal basis functions. The most 
common choice here has been the oscillator or Sturmian basis both of which 
have analytic Fourier-Bessel transforms making them well suited in calcula- 
tions where it is advantageous to work in configuration and momentum space 
simultaneously. In a variational Ritz-type approach the upper bounds of the 
true eigenvalues could be computed by diagonalizing the corresponding Hamil- 
tonian matrix (cf. [12], [13], [14]). The expectation values of the energy can be 
evaluated in momentum space and the potential expectation values in con- 
figuration space. The expansion method could be used in a similar fashion 
to solve the momentum space integral equation by means of the Galerkin 
method [15], [16]. With a judicious choice of the basis functions, the sin- 
gular integrals can be calculated analytically, or numerically. Note, that in 
this case the integrand is a known function and, therefore, the subtraction 
technique, like the one outlined in (2), is fully applicable. There are also non- 
variational approaches based on eigenfuction expansion such the collocation 
method [15], [17], or the Multhopp [4] [18] technique. Keeping N terms of the 
truncated expansion, the N expansion coefficients can be determined from the 
requirement that the integral equation be exactly satisfied at N distinct values 
of the momentum variable. The semi-spectral Chebyshev method developed 
in this paper also belongs to the last group. However, the Chebyshev series, 
after reshuffling takes the form of an interpolative formula. In consequence, 
the expansion coefficients and the function values taken at the mesh-points 
are connected by a linear relation (cf. [1]). Thus, put in a nut-shell, the un- 
derlying idea is to solve the integral equation exactly on the Chebyshev mesh 
and, subsequently, interpolate by means of a high degree polynomial. The 
plan of the presentation is as follows. In the next section we set the neces- 
sary background deriving the hypersigular integral equation associated with 
the Coulomb-plus-linear potential in momentum space. Upon introducing the 
Chebyshev mesh and using the interpolative formula for the wave function, 
the integral equation is converted into an algebraic eigenvalue problem. This 
is the ultimate form because the eigenvalue problem can be solved with the 
aid of standard library procedures. Section 3 is devoted to a numerical test 
where we compare the momentum space calculations with the results obtained 
by solving the Schrodinger equation in configuration space. Finally in the last 
section we present our conclusions. 
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2 Solution of the singular integral equation 

The Coulomb-plus-linear potential considered in this paper is V(r) = V^ c \r) + 
V iL \r) with 

V^ c \r) = -a/r; V {L \r) = r/a 2 (3) 

where the "coupling" a is dimensionless and the parameter a has a dimen- 
sion of length (h — c — 1 units are adopted hereafter). Both parameters are 
assumed to be provided. In momentum space the wave function <f>g(k) with 
orbital momentum £ obeys the partial wave Schrodinger equation 

(E - k 2 /2fx) Mk) = / Ve(k, k') (f>e(k') k' 2 dk' (4) 
Jo 

where fx is the quark- ant iquark reduced mass, E is the binding energy and 
Vg(k, k') denotes the £-th partial wave projection of the local potential V(r) 

2 r°° 

V e (k', k) — - j e (k'r) V{r) Je (kr) r 2 dr, (5) 

7T JO 

where je(x) is the spherical Bessel function [19]. Strictly speaking, upon in- 
serting (3) in (5), we obtain a divergent integral but a customary regularizing 
procedure to overcome this difficulty is first to multiply V(r) by a screen- 
ing factor e~ vr enforcing convergence and then set i] — > in the result. 
Applying this procedure, the Fourier transform (5) of a power-law potential 
v(r) = r 2n_1 , n — 0, 1, 2, ... can be effected in an analytic form [18] 

o r°° (2nV 
lim - / jAk'r) e-" r r 2n+1 jAkr) dr = , „ ~, Q n Az) (6) 

where z = (k 2 + k' 2 )/2kk' and the Qgiz) denotes n-th derivative of the Leg- 
endre function of the second kind with respect to the argument z (formula 
(5) in [18] contains a misprint). Setting n = and n — 1 in (6) we obtain, 
respectively, the kernels for the Coulomb (C) and the linear potential (L) 

V e {C \k,k') = -aQ e (z)/(7tkk'y, V e {L \k,k') = Q' e (z)/[7t(akk') 2 ]. (7) 

The Coulomb part of the kernel exhibits a logarithmic singularity for k' — k 
contained in the Legendre function. Indeed, the latter can be written as 

Q e (z)=P e (z)Qo(z)-w e - 1 (z) (8) 

where 

Q (z) = | log |(1 + z)/(l - z)\ = log \(k + k')/(k - k')\ (9) 

with Pg(z) being a Legendre polynomial. It is understood that the last term 
in (8) should be absent for £ = whereas for £ > it assumes the form of a 
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polynomial in z (cf. [19]) given by the expression 



we-i(z) = -Pn-i(z) Pe-n(z). 



n=l 



n 



(10) 



The kernel associated with the linear potential given in (7), in addition to the 
logarithmic singularity, exhibits also a second order pole, as may be seen by 
performing explicitly the differentiation in (8) 



with 



Q' e (z) = P[{z) Q (z) + P e (z) Q' (z) - w'^z) 



2kk' 



(11) 



;i2) 



l-z 2 \k' + k) (k'-k) 2 ' 
The second order pole given by (12) can be eliminated from the integral equa- 
tion (4) and to this end integration by parts is applied to this term. Quite 
generally, this procedure gives 



f 

Jo 



f(k,k') Mk') dfc' 
(k' - k) 2 



°° dk' d 



o k' — k dk 



- [f(k,k')Mk')] 



(13) 



where the unspecified function f(k, k 1 ) needs to be integrable. The above for- 
mula holds because the wave function <fii{k') vanishes when k' tends to either 
of the integration end points. The resulting Cauchy principal value integral 
in (13) can be computed by using the dedicated Chebyshev quadrature given 
in [1]. Nevertheless, the lowering of the order of the pole outlined above has 
its price and in the integral on the right hand side of (13) the derivative of 
the unknown wave function will appear. As we shall see in a moment, the 
semi-spectral Chebyshev method is well suited to handle such situation. 

It will be convenient for us using 1/a as the unit of energy, passing to dimen- 
sionless quantities: e = Ea, x = ka, x' = k'a. The resulting integral equation 



x 

2/ia 



<f>t(x) 



1 r°° ( 
— P/(*)log 

TlX JO I 



x' + X 



X 1 



X 



Wo 



4 
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dx' 
x' — X 



Xe(x') + <j) e (x' 



-IV 





(x) dx' 



x' 2 PAz) 





x' + X 


\p,{z) log 




x' — X 



dx' J (x' + x) 2 
w e -!(z) \ 4>e(x')x' dx' 



(14) 



involves two dimensionless parameters: a and 2/ia. Prime on a function of z 
denotes in (14) the derivative with respect to the argument. The derivative 
of the wave function appearing in the integrand of the second integral in (14) 
has been regarded as an additional function xe( x ) to be determined. In order 



6 



to complete our scheme the integral equation (14) needs to be supplemented 
with a complementary equation 

d<f>e(x)/dx = xe{x) (15) 

and we end up with two equations for two unknown functions: <f>g(x) and Xe( x )- 

The system (14)— (15) is amenable for computation and the integral equa- 
tion will be turned into a finite matrix equation. As a preliminary step, the 
semi-infinite domain of the independent variable x will be mapped onto a fi- 
nite interval (—1,1). Among endless possibilities perhaps the simplest is the 
rational mapping 

x = a(l + t)/(l-t), (16) 

where t G (—1,1) and a is a numeric parameter at our disposal providing 
additional control of the rate of convergence. We tried some other map- 
pings, specifically trigonometric (x = crtan[(7r/4)(l + £)]), or logarithmic 
(x = oTog[(3 + t)/(l — t)}) but they did not bring noticible improvement in 
the problem under consideration. The semi-spectral Chebyshev method uses 
Chebyshev polynomials as the basis functions. The Chebyshev polynomial of 
the first kind T N (t) of the order N is defined by the formula 

T N (t) = cos[A arccos(t)] (17) 

and has N zeros in the interval (—1,1), located at the points 

i i = cos[7r(i-i)/JV]; i = 1, 2, JV. (18) 

In the following the variable t will be discretized by using the classical Cheby- 
shev mesh (18) in which case N becomes the order of approximation to be 
selected by the user. The semi-spectral Chebyshev method interpolates the 
unknown function f(t) on the Chebyshev mesh (18) 

N 

/(0=E/(**)G f i(0, (19) 
1=1 

where Gi(t) denotes the cardinal function with the property G.i(tj) = 5ij. These 
functions can be constructed as superpositions of Chebyshev polynomials 

2 N 

G s (t) = -'£ , T i - 1 (t j )T i - 1 (t), (20) 

i=l 

where the primed sigma denotes a summation in which the first term should 
be halved. By taking advantage of the interpolative formula (19), the differ- 
entiation or integration of a function reduces to differentiation or integration 
of Chebyshev polynomials which in most cases is elementary and can be per- 
formed in an analytic form. In consequence, the array containing the values of 
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the derivative computed at the grid-points will be connected to similar array 
representing the function by a linear transformation 



l<m/U =T l D ij f(t j ), i = 1,2, ...,N (21) 

t=ti 3=1 



where is easily computed numerical matrix (cf. [1]). There are also various 
integration rules available. Assuming that the function f(t) is non-singular in 
the integration domain, we have 

1 N 

/ f(t)dt = Y / w t f(t l ), (22) 
J ~ 1 i=i 

which is Gauss-Chebyshev integration in which the weighting function is equal 
to unity. The weights Wi are all positive and their sum equals to 2. Similar rules 
can be derived for singular integrals. The Cauchy principal value integration 
can be performed using the automated quadrature rule 

r 1 m^ = £ Ui{T)mi (23) 

J-l I T i=1 

where it is assumed that r G< —1,1 >. When r coincides with either of 
the integration end-points the integral is undefined. The dedicated weighting 
functions ^(r) can be calculated analytically and exhibit logarithmic end- 
point singularity for r = ±1. Similar rule can be obtained for a weakly singular 
integral 

! N 

/ f(t)\og\t-r\dt = ^(r)f(U), (24) 

J ~ l i= i 
where it is assumed that r e (—1, 1). In contrast with the previous case, 
log \t — r\ singularity is integrable and the dedicated weighting functions fij(r) 
do exist even when r coincides with either of the integration end-points. For 
explicit analytic expressions for all of the weighting functions introduced above 
the reader is referred to [1]. 

To arrive at the ultimate finite matrix eigenvalue problem, as the first step, we 
map both, the external (x), and the internal (x 1 ) variable onto the (—1, 1) in- 
terval with the aid of (16). Subsequently, the problem is discretized by putting 
the external variable on the Chebyshev mesh (18), at the same time replacing 
all integrations in (14) by summations, following the appropriate Chebyshev 
rules listed above. In practice this procedure leads to a chain of substitutions 
to be made in the integrals occurring in (14), viz. 

x -> Xi = cr(l + ti)/ (1 - U); (j> e (x) -> (f>i(xi) = Xf, 

and 

x' -> Xj = cr(l + tj)/ (1 - tj); (f>e(x') -> (f>e(xj) = Xy, 
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where Xi are the unknown mesh values of the wave function to be determined. 
The derivative Xe( x j) * s eliminated in favor of Xj with the aid of the 
matrix, accounting for the change of variables 

(I -t) 2 N 

k=l 

Further substitutions associated with integration, respectively, are 

dx' -> 2awj/(l -tj) 2 , 



for non-singular integrals 



dx' 1 — ti 

— — p 



x' — x 



for principal value integral, and 

x' + x 



log 



x 1 — X 



dx' -> 2a 



Wj log 1 1 — ti tj I — flj (U 



for integrals involving logarithmic singularity. Finally, all integrations will be 
effected by carrying out a summation over j. It is worth noting that the 
diagonal terms i — j are always finite and all singularities are under control. 

When the indicated above manipulations have been accomplished, we end up 
with a homogeneous system of iV algebraic equations in which the N unknowns 
are the mesh-point values of the wave function (Xj) and the Schrodinger 
equation takes the desired finite matrix form 

Vv + J^Sy-e^ X, = 0. (25) 

The non-symmetric matrix Vij represents here the potential and results from 
evaluating the integrals occurring on the right hand side of (14) (the explicit 
form of is rather lengthy and will not be quoted here). When the kinetic 
energy term is lumped together with into a single matrix, eq. (25) presents 
a standard algebraic eigenvalue problem. If need arises, the non-relativistic 
Schrodinger equation (25) can be easily converted to the relativistic form (1) 
in the center-of-mass frame by changing just the kinetic energy term 




x 2 / (2/ia) — > \jx 2 + (ami) 2 + \J x 2 + (am 2 ) 2 - a (mi + m 2 ). 

Our calculational scheme is now complete and for assigned values of I and two 
dimensionless parameters s = 1/2/ia and a specifying the strength of the two 
potentials in (3), we are in the position to determine numerically the value of 
the binding energy e(£, s, a). In the particular case £ = and a = the exact 
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result is known and the binding energy is e(0, s, 0) = — s 2 / 3 z v where z v with 
v — 1, 2, 3, ... denotes a zero of the Airy function Ai(z) (cf. [19]). 



3 Numerical test 

We start the numerical test with the Coulomb bound state problem leaving out 
the first two integrals on the right hand side of (14). The hydrogen-like bound 
state problem in momentum space has already been considered in [1] but to 
determine the bound states we solved the secular equation. It is therefore of 
interest to repeat the Coulomb bound-state calculation in which the energy 
spectrum is obtained by solving the algebraic eigenvalue problem (25). The 
latter procedure is much simpler as there is no need to solve a transcendental 
equation. In all our computations we were using the linear algebra package 
LAPACK [20] as our eigenvalue solver. The results for the Coulomb potential 
are displayed in Table 1. Since in this case the exact eigenenergies are known 
analytically we present the absolute value of the relative error on each level 
as a function of the mesh size N. The nodal quantum number n enumerates 
the the different bound states for a fixed £ with n = corresponding to 
the ground state. We wish to recall that with non-symmetric matrices the 
accuracy of the standard library procedures is believed to be not as good as 
in the case of symmetric matrices. Nevertheless, as seen from Table 1, the 
convergence rate is exponential and N = 80 is sufficient for securing machine 
accuracy. There are not very many methods available that would be capable 
of achieving such a high precision. For comparison, in the last raw (entries in 
parenthesis) we give the relative error corresponding to the traditional method 
using the subtraction scheme [8] in which case the resulting eigenvalue problem 
is symmetric. The advantage of the semi-spectral method is manifest. 

As our next test we take on the linear potential setting a = in (25) and 
putting for simplicity s = 1 in our computations. The resulting binding ener- 
gies e for different £ values are displayed in Table 2 using the same conventions 
as in Table 1. For £ = 0, as the exact values we take the zeros of the Airy 
function tabulated in [19]. For £ > the values marked as exact have been 
computed by solving the appropriate Schrodinger equation in configuration 
space. For this purpose we used the ingenious algorithm developed in [21]. 
The code from [21] has been revamped for obsolescent features and the orig- 
inal Runge-Kutta driver advancing the solution from x to x + h has been 
replaced by a more accurate driver based on Chebyshev approximation as 
described in [1]. After the above changes, the typical relative error in all con- 
sidered here cases was estimated to be of the order of 10~ n . As a cross-check, 
we succeeded in reproducing the exact results for £ = up to eleven significant 
digits. To obtain the entries in table 2 for each £ value we needed to solve the 
algebraic eigenvalue problem (25) and in nearly all considered here cases we 
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managed to get seven significant figures which is more than adequate in all 
practical applications. Our results have been obtained keeping quite moderate 
approximation order N = 100. Only the £ = case which was more stubborn 
forced us to go to larger N. It is apparent from table 2 that the solutions are 
very stable with respect to increasing iV albeit the rate of convergence is no 
longer exponential. In fact, it is quite slow when compared with the Coulomb 
case. Making such comparison, however, it has to be kept in mind that in 
the linear potential case we need to determine two unknown functions (wave 
function and its derivative) rather than one and therefore N should have been 
doubled if we wanted to keep the same number of points per function. Other 
than that, there is probably a good deal of cancellation across the pole and 
this might be responsible for some loss of accuracy. 

Finally, we are going to consider the case where both, the Coulomb and the 
linear potential are present. The quark- ant iquark potential has been adopted 
from a realistic study [22] of charmonium (cc) and bottomium (bb) V(r) = 
—a/r + f3r where we stick to the parameter values provided in [22], namely 

a = 0.50667, (3 = 0.1694 GeV 2 , m c — 1.37 GeV, m b — 4.79 GeV. (26) 

The results of our computations are presented in Table 3. The quarkonium 
masses M displayed there have been obtained from the expression M = 2m q + 
E where m q is the quark(antiquark) mass. To determine the binding energy 
E the appropriate non-relativistic Schrodinger equation was solved in both, 
the momentum and the configuration space. As seen from Table 3 there is 
excellent agreement between these two approaches. 



4 Summary and Conclusion 

The aim of this paper was to demonstrate the strength of the semi-spectral 
Chebyshev method in solving integral equations whose kernels exhibit singu- 
larities of the Cauchy or the logarithmic type. Such equations may be en- 
countered in quantum mechanics as has been exemplified by considering the 
Coulomb-plus-linear potential bound state problem in momentum space. The 
latter problem is considered in this work for illustrative purposes and therefore 
we have gone in some details. The semi-spectral Chebyshev method has many 
advantageous features. First, it is very easy to use since it is based on a polyno- 
mial interpolation where both, the mesh and the polynomials, can be readily 
obtained in an analytic form. Second, the programming is exceedingly simple 
because differentiation or integration of polynomials can be performed analyt- 
ically and on a mesh these operations take the form of matrix multiplications. 
The presented method is well suited to handle singular integral equations 
(with Cauchy or logarithmic singularities) because automatic quadratures are 
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provided for evaluating singular integrals. This allows for a quick and seam- 
less discretization and since the integrals involving singular kernels have finite 
diagonal elements the Nystrom method is still applicable. Ultimately, the in- 
tegral equation is converted into an algebraic eigenvalue problem which can 
be solved directly by standard library procedures. There is no need to solve 
a complicated transcendental equation. Third, the method is highly accurate. 
This is because the approximation is global basing on a polynomial of a high 
degree. The eigenvectors contain the wave function values on the mesh and 
can be used to calculate various expectation values. If this is not enough, once 
the integral equation has been solved exactly on the mesh, the solution at an 
arbitrary point may be immediately obtained by interpolation. In conclusion, 
with the aid of the semi-spectral Chebyshev method the solution of a singular 
integral equation becomes no more difficult than the solution of a Fredholm 
equation. 

References 

[I] A. Deloff, Ann. Phys. (N. Y.), in press 
[2] R. Omnes, Nuovo Cimento 8 (1958), 316. 

[3] S.S. Chandrasekarn, Radiative Transfer, Dover, New York, 1960. 

[4] A. Robinson, Wing Theory, Cambridge U.P. Cambridge, 1958; 
K. Karamachetti, Principles of Ideal Fluid Dynamics, New York, 1966. 

[5] N. Brambilla et al., Heavy quarkonium physics, CERN Yellow Report, CERN- 
2005-005, available from the archive as hep-ph/0412158 

[6] E. Eichten and K. Gottfried, Phys. Lett. B 66 (1977) 286; E. Eichten et al., 
Phys. Rev. D 17 (1978) 3090; E. Eichten et al., Phys. Rev. D 21 (1980) 203 

[7] D. Flamm, F. Schoberl and H. Uematsu, Phys. Rev. D 36 (1987); S. B. Elegba 
and M. A. Rashid, Phys. Rev. D38 (1988) 2911. 

[8] D. P. Heddle, Yong Rae Kwon and F. Tabakin, Comp. Phys. Commun. 38 
(1985) 71; Yong Rae Kwon and F. Tabakin, Phys. Rev. C18 (1978) 932 

[9] A. Tang and J. W. Norbury, Phys. Rev. E63 (2001) 06670 

[10] T. W. Chiu, J. Phys. A 19 (1986) 2537 

[II] D. Eyre and P. Vary, Phys. Rev. D34 (1986) 3467 

[12] W. Lucha, H. Rupprecht, and F. F. Schooberl, Phys. Rev. D45 (1992) 1233; 
W. Lucha and F. F. Schoberl, Phys. Rev. A56, (1997) 139 

[13] L. P. Fulcher, Z. Chen and K. C. Yeong, Phys. Rev. D47, (1993) 4122; 
L. P. Fulcher, ibid. D50,(1994) 447 



12 



[14] L. J. Nickisch, L. Durand, and B. Durand, Phys. Rev. D30 (1984) 660; ibid D30 
(1984) 1995(E) L. Durand and A. Gara, J. Math. Phys. 31 (1990) 2237 

[15] J. R. Spence and P. Vary, Phys. Rev. D35 (1987) 2191 

[16] K. H. Maung, D. E. Kahana, J. W. Norbury, Phys. Rev. D47 (1993) 1183. 

[17] H. C. Jean, D. Robson and A. G. Wiliams, Phys. Rev D50 (1994) 5873 

[18] S. Boukraa and J. L. Basdevant, J. Math. Phys. 30 (1985) 1060 

[19] M. Abramowitz and I. Stegun, (Eds.), Handbook of Mathematical Functions 
(Dover, New York 1972) 

[20] E. Anderson, Z. Bai, C. Bischof, C. S. Blackford, J. Demmel, J. Dongarra, J. Du 
Croz, A. Greenbaum, S. Hammarling, A. McKenney and D. Sorensen, LAPACK 
Users 7 Guide, Third Edition, Society for Industrial and Applied Mathematics, 
(1999), Philadelphia, PA 

[21] P. Falkensteiner, H. Grosse, F. Schoberl and P. Hertel, Comp. Phys. Commun. 
34 (1985) 287 

[22] C, Quigg and J. L. Rosner, Phys. Rep. 56 (1979) 167 



13 



Table 1 

Relative errors on the computed Coulomb binding energies. The corresponding er- 
rors appropriate to traditional method based on subtraction are given in parenthesis. 



1 = 



N 


n = 




n = 1 




n = 2 


n = 3 


n = 4 


40 


4 x 10" 


12 


3 x 10- 


10 


4 x 10- 9 


3 x 10- 8 


2 x 10- 7 


60 


2 x 10" 


13 


1 x 10- 


n 


2 x 10- 10 


1 x 10- 9 


6 x 10- 9 


80 


2 x 10- 


14 


1 x 10- 


12 


2 x 10- 11 


1 x 10- 10 


6 x 10- 10 


80 


(4 x 10 


" 5 ) 


(9 x 10 




(2 x 10- 4 ) 


(4 x 10~ 4 ) 


(6 x 10- 4 ) 












e= i 






N 


n = 




n = 1 




n = 2 


n = 3 


n = 4 


40 


8 x 10- 


13 


2 x 10- 


13 


3 x 10- 9 


6 x 10- 8 


5 x 10- 7 


60 


2 x 10- 


14 


4 x 10- 


13 


3 x 10- 12 


1 x 10- 11 


2 x 10- 10 


80 


2 x 10- 


15 


4 x 10- 


14 


3 x 10- 13 


1 x 10- 12 


2 x 10- 12 


80 


(7 x 10 


-6) 


(5 x 10 


" 5 ) 


(2 x 10~ 4 ) 


(5 x 10- 4 ) 


(1 x 10- 3 ) 


N 


n = 




n = 1 




£ = 2 
n = 2 


n = 3 


n = 4 


40 


3 x 10" 




2 x 10- 


1U 


2 x 10- y 


5 x 10- 6 


8 x 10- 5 


60 


3 x 10- 


15 


6 x 10- 


1 A 

14 


1 x 10" 12 


7 x 10- 10 


2 x 10- 8 


80 


2 x 10" 


10 


2 x 10- 


lo 


4 x 10- 14 


4 x 10" 13 


6 x 10" 12 


an 


^1 X 1U 




(a v in 

l^D X 1U 


) 


Cx ^ i n~ 4\ 




(o 1 n — 3^ 

^ X 1U J 


N 


n = 




n = 1 




e = 3 

n = 2 


n = 3 


n = 4 


40 


2 x 10- 


9 


2 x 10- 


7 


9 x 10- 6 


3 x 10- 4 


3 x 10-' 3 


60 


1 x 10- 


12 


3 x 10- 


12 


8 x 10- 11 


2 x 10- 8 


6 x 10- 7 


80 


1 x 10- 


13 


2 x 10- 


12 


6 x 10- 12 


1 x 10- 10 


5 x 10- 10 


80 


(8 x 10 


" 5 ) 


(3 x 10 


" 5 ) 


(3 x 10~ 4 ) 


(1 x 1Q- 3 ) 


(1 x IO- 3 ) 
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Table 2 

Binding energy e for a linear potential. 

£ = 



N 


n = 


n = 1 


n = 2 


n = 3 


n = 4 


50 


2.338034 


4.087928 


5.520416 


6.786654 


7.943940 


100 


2.338099 


4.087947 


5.520543 


6.786702 


7.944111 


150 


2.338105 


4.087949 


5.520555 


6.786706 


7.944127 


200 


2.338106 


4.087949 


5.520558 


6.786707 


7.944131 


250 


2.338107 


4.087949 


5.520559 


6.786708 


7.944132 


300 


2.338107 


4.087949 


5.520559 


6.786708 


7.944133 


exact 


2.338107 


4.087949 


5.520560 


6.786708 


7.944134 








e = i 






N 


n = 


n = 1 


n = 2 


n = 3 


n = 4 


50 


3.361254 


4.884452 


6.207617 


7.405649 


8.515212 


100 


3.361255 


4.884452 


6.207623 


7.405665 


8.515234 


exact 


3.361254 


4.884452 


6.207623 


7.405665 


8.515234 








1 = 2 






N 


n = 


n = 1 


n = 2 


n = 3 


n = 4 


50 


4.248183 


5.629693 


6.868774 


8.009828 


9.075383 


100 


4.248182 


5.629708 


6.868883 


8.009703 


9.077003 


exact 


4.248182 


5.629708 


6.868883 


8.009703 


9.077003 








£ = 3 






N 


n = 


n = 1 


n = 2 


n = 3 


n = 4 


50 


5.050918 


6.331874 


7.504206 


8.593338 


9.632163 


80 


5.050926 


6.332115 


7.504646 


8.597127 


9.627263 


exact 


5.050926 


6.332115 


7.504646 


8.597117 


9.627267 
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Table 3 

Charmonium (cc) and bottomium (bb) masses (all entries in GeV) computed from 
the Coulomb-plus-linear potential [22] V(r) = —a/r + fir with the parameters given 
in (26). The upper (lower) values result from a calculation conducted in momentum 
(configuration) space using non-relativistic Schrodinger equation. In all momentum 
space computations the mesh size was N = 80. 



m c = 1.37 
n = n = l n = 2 



£ = 





3.0869 


3.6748 


4.1094 






3.0869 


3.6748 


4.1093 


i = 


1 


3.4988 


3.9544 


4.3388 






3.4987 


3.9543 


4.3388 


I = 


2 


3.7868 


4.1868 


4.5407 






3.7868 


4.1868 


4.5407 



m b = 4.79 
n = n = l n = 2 



I = 





9.4550 


10.0105 


10.3423 






9.4547 


10.0104 


10.3422 


I = 


1 


9.9171 


10.2582 


10.5318 






9.9170 


10.2581 


10.5318 


I = 


2 


10.1555 


10.4385 


10.6838 






10.1554 


10.4385 


10.6410 
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